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ABSTRACT 


With the advent of long term low Earth-orbiting satellites 
comes the requirement to maintain a specified orbital radial 
band. Optimal control theory suggests that periodic thrusting 
me. MOre efficient than forced Keplerian motion in obital 
maintenance through the use of Primer vectoring. This thesis 
examined the efficiency of fixed-angle transverse thrusting as 
an alternative to Primer vectoring. Numerical analysis shows 
that a thrust angle of 65-70 degrees is feasible for radial 
band control for a wide range of parameters. Fuel usage can 
be minimized through the proper selection of radial bandwidth, 
thruster size, and thrust angle. This thesis shows that 
forced Keplerian motion is always superior to fixed-angle 
transverse thrusting from a fuel usage standpoint, and hence 
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Le INTRODUCTION 


Orbital maintenance for low-Earth orbiting (LEO) 
satellites can be accomplished using forced Keplerian motion 
or by bang-bang control. Forced Keplerian motion utilizes 
constant thrust, with magnitude equal to drag, to prevent 
orbital decay and maintain orbital altitude. Bang-bang 
control allows the satellite orbit to decay and then fires 
thrusters to reboost the satellite and keep it within a 
prescribed radial band. This can be modeled as a series of 
orbital transfers offsetting decay due to atmospheric drag. 
Currently, LEO satellites are not required to maintain a 
specified orbital band, and are reboosted only to prevent 
reentry or to change orbit. 

Historically, co-planar orbital transfers have been 
optimally accomplished using two or three-impulse Hohmann 
maneuvers [Ref. l:pp 78-88]. These maneuvers are designed 
for one-time transfers between exo-atmospheric orbits, not 
for a series of maneuvers to offset drag and maintain an 
orbital band. Bang-bang control utilizes one thruster 
firing per reboost of the satellite, vice two or three 
firings used in a Hohmann transfer. 

In Ref. 2, the problem of minimum-fuel orbital 
maintenance for low-altitude, non-lifting bodies was 


considered, with the conclusion that fixed-angle, transverse 


thrusting (1.e. bang-bang control with the thrust veceoumem 
a constant angle relative to the local horizontal) is less 
etficient than forced Keplerian motion. This topic was 
first explored by Ross and Melton [Ref. 3] who concluded 
that forced Keplerian motion is not the optimal solution for 
this problem. Development of long-term, low-altitude 
satellites, such as the space station, dictate continued 
research into this topic to resolve this apparent 
COnlehod ler monn 

Little research has been done on the optimization of 
bang~bang control in relation to atmospheric effects on non- 
lifting (blunt) bodies. A better understanding of these 
effects is essential for optimization of orbits and 
propulsion equipment for long-term satellites, including the 
Space station. Additionally, design and operation of 


lifting bodies may benefit from this study. 


It. GENERAL FORMULATION OF THE PROBLEM 


MeO MeENOAS GS: OCDGiImizing Grbrecal transrers for LEO 
Satellites have been studied in Ref.2 and Ret. 3. On the 
basis of numerical simulations, Pauls [Ref. 2] concludes 
that bang-bang control is less efficient than forced 
Keplerian motion, and that transverse thrusting at a fixed 
angle of 70 degrees is generally the optimal feasible 
solution for orbital maintenance if bang-bang control is 
weed. his apparently contradicts the analytical results of 
Ross and Melton [{Ref. 3], who conclude that forced Keplerian 
motion is not optimal and that a bang-bang solution exists 
that optimizes fuel usage. This leads Pauls to conclude 
that variable thrust-vectoring is essential for optimality 
during finite-burns. However, Pauls study was limited to a 
small class of satellites (essentially space platforms) in 
the upper atmosphere. Therefore, his conclusions are not 
Globally valid, and questions remain. Is transverse 
thrusting less efficient than forced Keplerian motion for 
all parameters? What exactly are these parameters? Is the 
70 degree thrust angle always optimal? This thesis will try 
to answer some of these questions. 

The problem is defined as maintaining a spacecraft 


within a prescribed radial band R,,<R<R,,,- Comparison of 


fuel usage for bang-bang control and torced Keplerian 
motion, over a wide range of parameters, will provide 
insights on the optimality question. Also, determination of 
a thrust angle, or range of angles, that optimizes the bang- 
bang control solution is desired. 

Ross and Melton [Ref. 3] base their conclusions on 
methods utilizing optimal control theory, which require 
initial and boundary conditions to facilitate a solution 
uSing numerical integration. This is a cumbersome solution 
technique that requires recomputation for slight changes in 
initial or boundary conditions. While this may be the most 
accurate solution technique, it 1S not easy or very 
practical. 

Pauls [Ref. 2] develops a computer model to solve this 
problem. Equations of motion were derived and solved for 
various parameter combinations to optimize bang-bang 
control. Attempts to define a radial band and then control 
the spacecraft within that band were unsuccessful. However, 
maintaining specific energy allowed a band to be maintained, 
albeit larger than the desired band. All solutions to this 
model predicted that forced Keplerian motion is more 
efficient than bang-bang control. 

Further study using a computer model is dictated. Non- 
dimensionalization of the equations of motion will allow 
effects of variations in parameters to be more easily 


analyzed. Variations of parameters will be analyzed for 


trends and erfects on solutions, leading to optimal 
solutions for design of bang-bang control systems (sucn as 
PImMoter SlZe, wenrustmangle, and Size of radial band). 
Optimality will be based on fuel usage, maintenance of the 
required radial band (with minimai eccentricity), and 
Pm@G@usecer Liring patterns. An ability to predict the actual 
radial band controlled is desired. 

Results of the computer model will then be compared to 
those obtained using forced Keplerian motion to determine 
the optimal means of maintaining a spacecraft within a 


design radial band while subjected to continuous atmospheric 


olatste ie 


III. FORMULA DEVELOPMENT 


A. DEVELOPMENT OF THE EQUATIONS OF MOTION 

For simplification, orbital motion is assumed to be 
coplanar and initially circular. Assuming a non-1i1teme 
(blunt) body, drag becomes the only aerodynamic force 
perturbing the orbit and thrust is the external force 
applied to counter-balance drag. The external forces acting 
on the spacecraft consist of gravity, aerodynamic drags 
thrust. Figure 1 shows the coordinate system and net forces 


fomPehl st two=body prob Lem. 


Transverse axis 


Radial axis 





Figure 1 Graphical Representation of Coordinate System 


TRewequacilens OL MOtion Ler this two-body problem are 


a, =p? (2) 





where a. and a,, are the radial and transverse accelerations 
and © F. and Y F, are the sums of the external radial and 
transverse forces, respectively, and m is spacecraft mass. 


From Figure 1, the components of drag are 


D.=-Dsin(y) (3) 


D,=-Dcos(y) (4) 


= 
¢ 
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and the components of thrust are 


T =Tsin(a@) (5) 


Te eeec( a) (6) 


The angle y is the flight path angle and is defined by the 
intersection of the velocity vector and transverse axis, 
while angle a is the thrust angle and is defined by the 
intersection of the thrust vector and transverse axis. The 


equations of motion can now be written 


Oy ue ar C7) 
r? Mm Mm 

Bee pe Se entae e (8) 
m =m 


where uw is the Earth’s gravitational constant, @ is the 
angular position of the spacecraft, and primes denote 


differentiation with wespece, concime- 


B. NONDIMENSIONALIZATION OF THE EQUATIONS OF MOTION 

The equations of motion are nondimensionalized so that 
variations in parameters are minimalized. The 
nondimensionalizing constants are chosen to "balance" the 
equations and thus produce the same order of magnitude 
responses. This facilitates study of the effects of 
changing parameters on the optimization problem. 

1. Definitions 

In order to nondimensionalize the equations, 

variables must be defined and nondimensionalizing constants 
must be chosen. Nondimensionalizing constants for radius 
and time are designated R and T, respectively, and were 
chosen somewhat arbitrarily based on the Earth’s radius. 
Thus, R is the radius of the Earth and Tt is the period of a 


circular orbit whose radius equals R, and is defined as 





r=an| & 2) 
U 


These constants define the nondimensionalized variables for 


radius and time 


ty| 
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The nondimensionalizing factor for mass is M (spacecraft 


initial mass), wnich defines the nondimensionalized mass 


a= (25 
M 
The variable #@ is already dimensionless. 
2. Nondimensionalization 
Using the above definitions, the following 
relationships are developed. 
art ae (13) 
bE » Us 18RE 
= 2 42 
epee oe (14) 
Mie Rk ae R qt? 
dé _, a (15) 
dt ae 
2 2 
add (7 WP) od (16) 
die aa “Ge aie 


Substituting these equations into the equations 7 and 8, 


substituting equation 9 for t where convenient, and 


rearranging terms, yields the following nondimensionalized 


results 











=r(—_)- ane ae (ize 
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lee Acie ie Wiseey , Bo Meese (18) 
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C. NONDIMENSIONALIZATION OF ATMOSPHERIC DRAG 

1. Definitions 
Drag is given by 
D=+ps_cw? (19) 
xP rep—d 


where p is atmospheric density (altitude dependent), S,, is 
the reference surface area of the spacecraft affected by 
atmospheric density, C, is the coefficient of drag, and Vv is 
spacecraft velocity. 
2. Nondimensionalization 
Nondimensionalization of v’ is accomplished using 


the relationship 


Sb as dé 
2 Mee fh) ae Hy 2 capa) z (20) 
veer sro (ae) Nee! 


Substituting equations 10,13, and 15 gives 


J8, 


2 rm, 5.4 | Zon 
ee ree (21) 
ia qd Ce 4 tat 





which is then substituted into equation 19 and simplified 


into the form used in equations 17 and 18 to give 


aap ee PS CGR V _ pRV- (22) 


RmM R  mMr2 2mB 


mere 8B 15 the initial ballistic coefficient and defined as 





B= : C23) 


This result is then put back into equations 17 and 18 to 


yleld the nondimensionalized equations 


p . 
aioe -F( ames _pRvsiny, T ane le (24) 
de 2mB RM 
2 2 2 
Gabe dédr_pRV’cosy, T° Tcosa (25) 





de £F dtdt 2MBr RE iM 


Inspection of these equations reveals that further 
Simplification is possible. Assuming an exponential 
atmospheric density model (which will be used throughout the 
study of this problem), with atmospheric scale height £6 


defines 


Dae eae (26) 


iis allows the ballistic coefficient, B, to be 


nondimensionalized as 


al 


B 
Oi 


is 





(27) 


The thrust terms can be similarly nondimensionalized by 


defining a nondimensionalized thrust 


(28) 





Substituting equations 27 and 28 into equations 24 and 25, 
and nondimensionalizing the exponent for drag result in the 


final equations 


-BR(E- Tee) 
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These equations will be used to generate the computer model 
to study the optimality problem of bang-bang control. The 
only parameters that need to be varied to study the 


optimality problem are 


On ae 


Since the only variable parameter in the nondimensionalized 
thrust term is thrust, T, either one can be varied for the 


study. 
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D. DEVELOPMENT OF MASS EQUATIONS 


One or the optimality criteria to be used in this study 





Uemeene mass of fuel consumed to maintain orbit. The 
melactivonsnip 
dm: re 
Oem ~ sp 


defines the rate of fuel consumption and can be integrated 
over time to determine the total fuel consumed. Using 
equations 11, 12, and 29, this can be nondimensionalized and 
Simplified to 


dm __ R 7 
cee 2 eG 





(32) 


These equations will be used to calculate the fuel used to 
maintain orbit using bang-bang control for comparison to 


fuel used in forced Keplerian motion. 


ies 


IV. DEVELOPMENT AND VALIDATION OF COMPUTER MODEL 


A. COMPUTER PROGRAM DEVELOPMENT 

A computer program, written in FORTRAN and listed in 
Appendix A with sample input and output files, was developed 
to simulate spacecraft orbital motion. A fourth-order 
Runge-Kutta numerical integration routine is used to 
integrate the equations of motion. The program is comprised 
of a main program section, which controls input, outpueameme 
flow of information, and five subroutines that provide 
computations necessary to Simulate orbital motion and 
Calculate the osculating orbital parameters. 

The non-dimensionalized equations previously developed 
(equations 29 and 30) are simplified for use in the computer 


program using the following state variable definitions. 


X,=r (33) 
= = le (34) 
Waites ae 

x,=6 (35) 
-% _ ad (36) 
eae ate 


These state variables are nondimensionalized as follows 
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ae (27) 
jy OM 3. 
ee (38) 
“chee ae 
y= WOW, GO ry (32) 
dé at 
The state variable x, is already nondimensionalized. These 


equations are then substituted into equations 29 and 30 to 


Give the equations of motion used in the computer program. 
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The first subroutine, DRAG, calculates the atmospheric 
drag experienced by the spacecraft. A constant atmospheric 
density model is used for initial program validation, after 
which an exponential density model is used. The second 
subroutine, EQN, updates the equations of motion that define 
the spacecraft’s orbital motion, equations 40 and 41. The 
third subroutine, ORBPAR, calculates the osculating orbital 
parameters for the spacecraft’s motion, including semi-major 
axils, apogee, perigee, and period. The fourth subroutine, 


RK4, is a standard fourth-order Runge-Kutta numerical 
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integration routine that integrates the state variables. 
The last subroutine, THRUST, governs ehe activation ana 
deactivation of the thrusters to maintain the satellite’s 


SE Dieta. 


B. PROGRAM VALIDATION 
1. Initial validation with no external forces 

Program validation proceeded in a logical series of 
steps, with each step requiring positive validation before 
proceeding to the next step. Initially, all external forces 
except gravity were neglected in order to maintain a 
circular orbit. As shown in Appendix B, Figures 1-4, 
radius, velocity, angular momentum, and specific energy all 
remained constant. Initial conditions of an elliptic orbit 
were then input into the program with all external forces 
except gravity neglected. Figures 5-8 of Appendix B show 
that semi-major axis (SMA), eccentricity, specific energy 
and angular momentum are constant. This validates the 
Initial) programe 

2. Validation of drag 

A constant atmospheric density was then introduced 
into the program, with a value that would produce a 
noticeable effect in a small number of orbits to enhance the 
validation process. The program was run for ten orbits and 
the results compared with analytical results obtained using 


equations from Ref. 4, 
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oan. 
Aa=-2m( ——) pa- (42) 
See. 
Av=-1(—") pav (43) 
ee ae - 
eee (ee (44) 
ru V 





which show the changes in semi-major axis, velocity, and 
Perevoa per orbit of a spacecraft under drag. In all of 
these equations, ballistic coefficient (B) was substituted 
po M/CjS... The program was run Using both radius and semi- 
major axis for the term "a" in the above equations. Results 
should be very similar since semi-major axis and radius are 
assumed to be equal in the formulation of the analytic 
equations and should be relatively equal for the small 
changes expected. Percent errors were calculated for the 
computed vs analytical results for changes in semi-major 
axis (SMA), velocity (VEL), and period (PER), as well as for 
the percent differences between computed and analytical 
values of semi-major axis, velocity, and period after each 
orbit. Table I shows the results when radius is used in the 
analytical equations for "a", while Table II shows the 
results when semi-major axis is used in the analytical 
equations for "a". As shown in Table I, the percent error 
for the changes in semi-major axis (%ASMA error) and period 


(sQPER error) start small after one orbit, grow through the 


iby, 


Table I COMPUTED AND ANALYTIC ERRORS USING RADIUS 






















ORBIT | <ASMA SMA SAP % PER SAV % VEL 
7 ERROR — ERROR ERROE ERROR ERROR 


1 90000 | 7.2875 | .00015 
1.6709 | .00003 | 1.6710 | .00001 | 5.1534 | .00029 
00043 
| 6 | 10.037 | .00007 | 20.027 | .o0022 | 10.484 | .o0081_ 
11.860 | .00010 | 11.860 | .00018 | 14.265 | .00107 | 
SEAR. 13.859 | .00120 




















Table II COMPUTED AND ANALYTIC ERRORS USING SEMI-MAJOR AXIS 











oo S6ASMA % 6 AP % PER 6AV VEL 
eae nati ERROR ERROR ERROR ERROR 


Be oe ) 9.5531 -Q00001 | 9455531 | .00002 | 722875 Asc CG 
3 






2 00029 
[3 [9.9531 | .o000a | 9.5532 | 00006 [1.0289 | .o0043 
[5 [9.5532 | .00006 | 9.5534 | .00009 | 6.6409 | .00069_ 
| 6 | 9.5532 | .00007 | 9.5535 | .o0021 | 10.484 | 00082 | 
9.5533 | .00010 | 9.5536 | .00015 | 14.265 | .00107 
C2 Tese nfs ann [son 














18 


cifth orbit, and then steady out at roughly eleven percent. 
This final error vercentage closely corresponds with the 
relatively constant 9.5% error found when using semi-major 
mec in Table IT. However, in both cases, aS shown in 
fmieres Lt and I7, the percent error between the calculated 
and analytic semi-major axls (% SMA error) and period (% PER 
error) after each orbit is less than 0.0003%. The absolute 
percent error in the change in velocity (#AV error) and the 
difference between the calculated and analytic velocity (% 
VEL error) after each orbit are the same whether using seml- 
major axis or radius for "a", as shown in Tables I and II. 
Again, while the percent error in the change in velocity 
Grows to greater than ten percent, the percent error between 
calculated and analytic velocity after each orbit very small 
(less than 0.002%). The reason for the ten percent errors 
between parameter changes appears to be the accuracy 
involved when the changes are so small (i.e. semi-major axis 
1s changing by approximately 200 meters per orbit out of an 
initial orbit of 6638200 meters). Thus an error of 10% in 
change of semi-major axis is only 20 meters, which is 
insignificant in terms of the initial radius. The analytic 
and computed results are thus in agreement and the drag 


portion of the program is validated. 


ro 


3. Validation of thrust 
The thrust portion of the program was validated 
throughout the study by having drag present and the 
Spacecraft maintain an orbital band by reboosting using 


EnEFUSters. 


Ge ORBIT CONTROL STRATEGY 

It is desired to design a control strategy that is 
Simple and will maintain the spacecraft within a pre- 
specified radial band. Figure 2 shows a typical band, 


defining limits of radius and specific energy. 





Figure 2 Orbital band with radius and energy 
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All control strategies were run using numerous combinations 
of parameter values (i.e thrust vector angle, a, was varied 
from 0-90° in 5° increments, radial band was varied using 2, 
Zomana 100 km bands, thrust was varied from 25 to 300 N., 
etc.). For consistency, all cases presented here were run 
with the following baseline values 

RO=6638.2km 

M=20000kq 

Thrust=300N 

Ballisticcoefficient=B=150kgq/m- 


9,=9.4e-10kg/m°? 


Gi aepeeie =e Kem 
2 sp-300sec 


“%=70° 
Radialband=2km 


which produce representative results corresponding to 


B=25000 
T=0.060 


1. Control using radius 

Initial attempts were made to control the radial 
band using radius as a parameter, i.e. turn thrusters on 
when the radius drops to the low end of the band (RMIN) and 
turn them off when the satellite is reboosted to the high 
end of the band (RMAX). As shown in Figures 1 of Appendix 
C, this control strategy was ineffective because the 
specific energy of the orbit continually increased and never 
leveled off (see Appendix C, Figure 2). The eccentricity of 


this orbit grew to approximately 0.1. 


ae 


2. Control using total energy 
Control was next attempted using total energy. The 
total energy of the initial orbit was calculated, as were 
the total energies of the spacecraft in circular orbits at 
the minimum and maximum of the radial band (see Figure 2) 
using 


gy=m( 2-H) (45) 


The thrusters were turned on when the spacecraft total 
energy dropped below the minimum orbit total energy and 
turned off when the spacecraft total energy rose above 
maximum orbit total energy. This method failed to control 
radius, as the spacecraft orbit always eventually decayed, 
but successfully controlled energy, as shown in Figures 3 
and 4 of Appendix C. While this method succeeded in 
maintaining a pre-determined radial band, the orbital radius 
was not maintained and decayed. 
3. Control using radius and total energy 

A combination control strategy utilizing both radius 
and total energy was attempted. MThrusters were turned on 
when the orbital radius dropped to the minimum radius and 
turned off when total energy increased to the initial total 
energy. Figures 5 and 6 of Appendix C show the unsuccessful 
results. Total energy would seem to be the parameter of 


choice for maintaining the orbit vice specific energy. 
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However, if the total energy 1S maintained at a constant 
value the orbital radius must decrease and/or velocity must 
decrease Since the spacecraft mass decreases due to 
thrusting. This can be shown using equation 13. 
4. Control using specific energy 

The next attempt at control involved spacecraft 
specific energy. This control was similar to the attempt 
uSing total energy, only specific energy was used. Specific 
energy 1S a function of orbital radius (for circular 
orbits), so "conserving" specific energy should result in a 
"conservation" of radius (l.e. a constant average radius). 
As shown in Figure 7 of Appendix C, the orbital radius was 
maintained, however the orbital band was not. Conserving 
Spacecraft specific energy, as shown in Figure 8 of Appendix 
C, failed to maintain the initial radius of the spacecraft. 

5. Control using radius and specific energy 

Control using radius and spacecraft specific energy 
is successful in maintaining radius and a radial band, as 
shown in Figures 9 and 10 of Appendix C. Thrusters are 
turned on based on the spacecraft orbit decaying below a 
specified radial minimum and spacecraft specific energy 
being below the initial specific energy. Thrusters remain 
on until specific energy increases above the initial value. 
This control strategy results in an orbit whose radius 


oscillates about the initial radius, but not necessarily at 
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the pre-specified value. The variation of the actual radial 
band from the desired band varies and iS examined in greater 


detail in the next chapter. 
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V. ANALYSIS AND RESULTS 


Analysis of the validated control strategy is performed 
in five steps. A set of parameter values is chosen as the 
baseline case about which they will be varied. These values 
are selected as representative of a space-station/platform 
type satellite. The baseline values are 

R=6638.2km 
M=20000kg 
Jep-300SeC 
T=300N 
9,=9.4e-10kg/m°? 
B=150kg/m°* 
B=25000 
ieee Sian * 
The simulation program is run for 100 orbits, with the 


equations of motion updated 5000 times per orbit, and output 


sampled 100 times per orbit. 


A. VARIATION OF THRUST ANGLE, @ 

First, numerous cases are studied to determine the 
validity of an optimal thrust angle, a, of 70°. The program 
is run with combinations of variations of parameters about 
the baseline values, while varying a from 0-90° in 5° 
increments. In all cases, thrust angles of 65 and 70° are 
the only values that maintain a constant radial band. 


Figures 1-4 of Appendix D show results over the range of 60- 
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75° for the baseline parameters, and are representative of 
all cases. Values of a@ below 60° and above 75° vroduce 
results that degrade significantly from those shown. Semi- 
major axis (SMA) and eccentricity are also compared 
throughout the analysis for determining optimal values of a. 
Figures 5 and 6 of Appendix D show results consistent with a 
maintained radial band (a=70°), while Figures 7 and 8 show 
results consistent with a radial band not being maintained 
(~@=50°). As can be seen, eccentricity must level off toa 
constant value to maintain a radial band, while SMA alone is 
not an indicator of orbital maintenance. 

The mass of fuel required to maintain orbit generally 
increases as a increases from 0 to 90°, with some minor 
local fluctuations. Table III shows the mass of fuel burned 
in kilograms (kg) to maintain 2, 25, and 100 km pre- 


specified bands, with @ varying from 60-75° and baseline 


Tanlte TIT MASS OF FUEL BURNEDE ns) 


Band (kn) 


2645 2013 2174 











values for other parameters, and is typical of the trend 


seen. Thus, from a fuel usage standpoint a should be 
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minimized, but from a radial band standpoint a must be in 
Pao waaige Of 65-70", 

This shows that a constant thrust angle, a, of 
approximately 65-70° is necessary to maintain a radial band. 
The following analyses are conducted with a range of 60-75° 
for a to minimize the data collection while ensuring the 


Continuing validity of this result. 


B. VARIATION OF PRE-SPECIFIED RADIAL BANDWIDTH 

Once a is chosen, it is desired to know how the pre- 
specified radial bandwidth affects the results. While 
maintaining @ within the 60-75° range and using the baseline 
values for other parameters, the program 1s run using 
Beeec tied Fadial bandwidths of 1,2,5,10,25,100, and 200 km. 
Although the large bandwidths are not practical, they are 
used to determine the existance of any limiting conditions. 
These results, shown in Table IV (first entry is fuel burned 
in kg, second entry is actual radial band maintained in km), 
show that fuel required to maintain the orbit generally 
increases as the controlled radial bandwidth increases, with 
the exception being the 100 km band. Although the error 
between the radial band maintained and the pre-specified 
band is minimized as a@ increases, the maintained band cannot 
be correlated to the pre-specified band (changing any 
parameter changes the maintained band). For this reason, 


the radial bands discussed throughout this analysis are pre- 


PAS 


Table IV MASS OF FUEL BURNED (kg) ,AND ACTUAL RADIAL BAND 
MAINTAINED (km) 





—— = =_ * 


Band a=60 a=65 





a=70 | a=75 
| 1 tS 3 oO ay HoS4 Ong 202.5719 o 262 /7 ee 
2 2049/2 20137 leno | 2174/1. 3 268 372 


2516/12 2628/8 2583/6 3065/5 
2392/25 2554/16 3191/13 3602/11 
Aas 2643/55 2828/43 eee 3595/30 


| 
2812/90 3049/75 | 3136/65 4289/58 


3020/130 | 2859/125 | 2826/60 3787/75 
200 2993/130 | 3549/140 | 4470/120 | 5743/135 
















specified bands. Figures 1-3 of Appendix E show the mass of 
fuel burned using bang-bang control and forced Keplerian 
control for 2, 25, and 100 km bands with a=70°. As can be 
seen, the fuel used increases in a generally linear manner 
as the number of orbits increases. For the 100 km band, the 
fuel used is linear when the points at the end of each burn 
are connected, as indicated by the dotted line in Figure 3. 
Thus, these results can be extrapolated for any number of 
orbits to give fuel usage vs Keplerian usage. However, the 
thrusters do not fire at the same point in each orbit (the 
thrusters fire approximately every three-fourths of an orbit 
for the cases of a small radial band and, thus firing 
precesses each orbit). These results are consistent 


throughout the analysis. 
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A constant radial band 1s maintained for a@ of 
approximately 65-70° but not for 60 and 75°, aS snown in 
Figures 4-7 of Appendix E, which again verifies the required 
thrust angle. These figures are for a pre-speciried radial 
bandwidth of 25 km, the results are similar for bandwidths 
Sie oO Km. Figure Seer Appendix E shows the results for a 
pre-specified bandwidth of 100 km and an @ of 70° (this is 
representative of results for a=60-75° and bands of 100 km 
and greater). As can be seen, the band is wide enough that 
PiewoLems fire infrequently to reboost the satellite. The 
actual radial bandwidth maintained was then compared to the 
specified band. Based on these results, all subsequent 
analysis is performed for three cases of radial bandwidth; 
ee2e>, and 100 km. This will give results for a small, 
medium and large bandwidths, while examining the unusual 


results at 100 km. 


C. VARIATION OF ATMOSPHERIC SCALE HEIGHT, £8 

Of the three variable parameters in the non- 
dimensionalized equations, ~@ should have the most 
predictable effect on the results and is analyzed first. 
Using baseline values for other parameters, 6 is varied from 
1.0e€-10 to 1.0e-2, and results are examined for logical 
trends. Logic would dictate that increasing the value of £6 
Will increase the rate of orbital decay within the band, 


requiring more thruster firings to reboost the satellite, 


Za 


and thus require more fuel. This is indeed the case as the 
results in Appendix F show. Figure 1 shows the fuel used as 
a function of $ for a pre-specified 2 km band, while Figures 
2 and 3 show the same results for 25 and 100 km bands, 
respectively. All three plots are similar in appearance 
indicating changes in $ affect orbital decay uniformly over 
the range of bands. Below a certain value, le-5 for 2 km 
band, changes in 6 have no effect on fuel usage. Above a 
certain value, le-3 for 2 km band, the orbit decays too 
quickly for thrusters to control. In between these values, 
increases in 8 increasingly affect the fuel used. Since 
changes in @ do not affect the forced Keplerian solution, 
the mass of fuel required for forced Keplerian maintenance 
for all cases is approximately 649 kg. Thus, smaller values 
of 8B allow the bang-bang control to be closer in efficiency 
to forced Keplerian motion. Analysis of the plots for 
radial band again show that a radial band is not maintained 
for 60 or 75° thrust angle, but is maintained and optimized 
in the 65-70° range. From these results, it is shown that 
changes in 6 affect the efficiency of bang-bang control as 


expected and are predictable. 


D. VARIATION OF THRUST, T 
The next variable to be analyzed is thrust (T). As 
explained in Chapter 3, either thrust or nondimensionalized 


thrust can be varied and the results will be the same. 
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Thrust is chosen because variations in thrust are more 
Pact, related tO and wnderstood. For this analysis, thrust 
1S varied over the range 10-1000 N, while all other 
parameters are fixed at the baseline value. Figures 1-3 of 
Appendix G show the fuel required, for a=70°, to maintain 2, 
25, and 100 km pre-specified bands, respectively. 

Trends are similar on all three graphs. Fuel usage 
Srertes Off low for small values of thrust, with thrust being 
small enough to produce quasi-forced Keplerian motion 
(thrusters fire for majority of time to reboost satellite). 
Forced Keplerian thrust is determined to be approximately 
3.8 N for the reference conditions. To obtain this value of 
transverse thrust using a thrust angle of 70° would require 
a thrust of approximately 11.1N. Thus, small values of 
thrust may produce a quasi-forced Keplerian motion with 
respect to the transverse thrust component. This is 
followed by a period of increased fuel usage as thrust 
increases. In the 2 km band case, increasing thrust further 
shows a drop in fuel usage followed by a relatively constant 
fuel usage as thrust increases above 300 N. The 25 km band 
case shows more fluctuation in fuel usage as thrust 
increases, followed again by a relatively constant fuel 
usage as thrust increases above 300 N. For the 100 km band 


case, increasing thruster size generally increases fuel 


ob 


usage except for a fairly constant region in the rangewem 
DOS 2 SOn Nic 

Thus, other than a possible minimizing of fuel usage due 
to this apparent quasi-forced Keplerian motion, changes in 
thrust appear to have little effect on the optimality 
problem. The smaller fuel usage totals at low thrust (25-50 
N) may be indicative of approaching a quasi-forced Keplerian 
motion (except thruster angle is 70° vice 0°). For all three 
cases, as thrust increases, fuel usage appears to approach a 
constant value. An optimal thrust angle of 65-70° is again 


verified in this portion of the analysis. 


E. VARIATIONS IN NONDIMENSIONALIZED BALLISTIC COEFFICIENT 

The last parameter to be analyzed is nondimensional 
ballistic coefficient, B. By nondimensionalizing, 
variations in B will encompass variations in the individual 
parameters of mass, atmospheric density, coefficient of 
drag, and reference surface area, as well as combinations of 
Variations involving these parameters. 

Analysis is conducted by varying B over the range 4e3- 
2e5 and comparing mass of fuel burned using bang-bang 
control (MFB) to mass of fuel burned using forced Keplerian 
motion (MFBK). This analysis is conducted for the three 
bands of 2, 25, and 100 km using baseline values for other 
parameters while varying a from 60-75°. Figures 1 and 2 of 


Appendix H show the results for the 2 km band case, results 


2 


for the 25 and 100 km bands are similar. As seen, the mass 
eee “el burned for both bang-bang control and forced 
Keplerian motion decreases as B increases, approaching a 
constant value as B increases above 2e5. FOE Optimalicy, 
maximizing B minimizes fuel usage. 

However, this does not give a feel for how erficient 
bang-bang control is compared to forced Keplerian motion. 
To analyze this, the ratio of MFB to MFBK is plotted versus 
B for 2, 25, and 100 km bands and @ of 65 and 70°. These 
results are shown in Figures 3-8 of Appendix H. Again, B is 
varied over the range 4e3-2e5. Values of B below and above 
this range (and the range plotted for 100 km band) cause the 
orbit to decay too rapidly or too slowly, respectively. 
Logically, this is analogous to increasing the atmospheric 
density to a value that causes the orbit to decay faster 
than thrusters can maintain, or decreasing the atmospheric 
density to a value that causes the orbit to decay so slowly 
that the thrusters never fire in 100 orbits (thus giving no 
data since at least one firing of the thrusters in bang-bang 
control is required for comparison to forced Keplerian 
mocion) . 

As seen in Figures 3-8, the general results for a of 65 
and 70° are similar for all three bands. Also, the results 
for all three bands show a roughly constant value for the 


ratio of MFB to MFBK, with the 2 km band being most 
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efficient compared to forced Keplerian motion, followed by 
the 100 km band and then the 25 km band. A value or ieee 
the ratio would indicate that the bang-bang Control wae 
equally efficient to Keplerian motion. Values less than 1 
would indicate that the bang-bang control 1S superior to 
Keplerian motion and would numerically validate the 
theoretical results of Ross and Melton [{Ref. 3]. The 
results shown are consistently above a ratio of 3, 


indicating that forced Keplerian motion is significantly 


eel 


more efficient than bang~-bang control over all values of 
Throughout this analysis, the optimality of a in the 

range of 65-70° was verified. This result has remained 

valid throughout all of the analyses conducted, indicating 


that this is a globally valid result. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 


The purpose or this thesis was to determine if there are 
any cases where fixed angle, bang-bang control, orbital 
maintenance approaches or exceeds the efficiency of forced 
Keplerian motion as proposed by Ross and Melton [Ref. 3]. 

By nondimensionalizing the equations, a thorough study of 
the effects of all parameter combinations was possible. 
Analysis of the results indicate that, in all cases, forced 
Keplerian motion is superior to bang-bang control. There 
mec nO Cases where bang—-bang centrol even remotely 
approaches the efficiency of forced Keplerian motion. Thus, 
these results appear to be globally valid. 

The optimality question for thrust angle also appears to 
be resolved globally. In all cases, the thrust angle is 
required to be between 65 and 70° in order to maintain a 
constant radial band. Generally, minimizing this angle 
improves fuel consumption, so the angle should be minimized 
as far as practicable while ensuring a constant radial band. 

Maintenance of a pre-specified radial bandwidth is not 
always possible or predictable for this control strategy. 
The error between specified and maintained band can be 
minimized through proper bandwidth selection. Selecting a 


relatively small, 2 km, or large, 100 km, band minimizes 


a0 


the error. Also, maximizing the thrust angle, a, generamm 
minimizes the error, and thus, a should be set as high as 
possible, while maintaining a constant radial band. Since 
this contradicts the requirement to minimize a@ for fuel 
savings, trade-offs must be considered. Minimizing the 
specified bandwidth moves the efficiency of bang-bang 
control closer to that of forced Keplerian motion. 

The effects of changes in atmospheric scale height, 6, 
are logical and predictable on the effects of the mass of 
fuel burned in bang~bang control. Since, by definition 
changes in @ do not affect forced Keplerian results, the 
changes in mass burned in bang-bang control directly 
correlate to changes in efficiency. Thus, the lower the 
value of 6, the more efficient bang-bang control is in 
relation to Keplerian motion. 

The effects of varying thrust are not so predictable; 
however, general trends can be inferred from the results. 
Low values of thrust appear to approach quasi-forced 
Keplerian motion and are more efficient. Higher values of 
thrust increase the mass of fuel burned, however, a maximum 
limit appears to be approached. As in the case of £8, 
changes in thrust do not affect the mass of fuel burned in 
Keplerian control, thus the mass of fuel burned in bang~bang 
control is indicative of the efficiency. 

This is not the case for changes in B, where changes 


affect the final Keplerian fuel usage. Thus, in order to 


36 


examine erficiency, the ratio of mass of fuel burnea using 
pana-bang control to mass of fuel burned with Keplerian 
motion must be used. It was initially expected that changes 
in B would have noticeable effects on efficiency and would 
show a trend that would indicate values where bangq-bang 
control erficiency approached Keplerian efficiency. This is 
Memetne case, as changes in B seem to have little effect on 
the ratio of efficiencies. 

Thus, the superiority of forced Keplerian motion to 
fixed angle bang-bang control is globally confirmed. 
However, further research should be conducted into the 
solutions proposed by Ross and Melton [Ref.3]. Control 
uSing variable thrust angle bang-bang control (so-called 
Primer-vectoring) should be examined for comparison to these 
results. Also, fixed angle bang-bang control, as described 
here, should be compared with other innovative methods of 
control, as well as with the "Lambert control" used with the 
shuttle, for possible fuel savings for space-station type 
e2actiorms. 

The results of this thesis also suggest areas for 
further study. A method for maintaining a pre-specified 
band is highly desirable, but was not achieved in this 
study. Although the global validity of the thrust angle 
being required to be between 65 and 70° was shown, the 


reason for this has not been studied. The problem of 


ce 


micro-gravity constraints and its effect on thruster firings 
and control strategies has also not been studied, but would 
be a critical factor in any control design. The facumemeas 
rorced Keplerian motion has been shown to be on the order of 
three times more efficient than bang-bang control indicates 
that there is much room for research into further 
Optimization techniques for bang-bang control. Finally, the 
question of practicality has not been addressed in relation 
to both forced Keplerian motion (generally not very 
practicable but very efficient) and bang-bang control (very 


practicable but not very efficient). 
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APPENDIX A 


PROGRAM LISTING 


OBJECTIVE: COMPARISON OF FORCED KEPLARIAN TO BANG-BANG 


CONTROL OF ORBITAL ALTITUDE 


VARIABLE DEFINITIONS 


XBAR(1)=NONDIMENSIONALIZED ORBITAL RADIUS 
XBAR(2)=NONDIMENSIONALIZED ORBITAL RADIAL VELOCITY 
XBAR(3)=THETA 

XBAR(4)=NONDIMENSIONALIZED ORBITAL ANGULAR VELOCITY 
XBDOT(1)=DERIVATIVE OF XBAR(1) 

XBDOT(2)=DERIVATIVE OF XBAR(2) 

XBDOT(3)=DERIVATIVE OF XBAR(3) 

XBDOT(4)=DERIVATIVE OF XBAR(4) 

RO=INITIAL ORBITAL RADIUS 

R=ORBITAL RADIUS 

RREF=REFERENCE ALTITUDE FOR DENSITY 

D=DRAG (N) 

DK=DRAG KEPLERIAN 

E=SPECIFIC ENERGY 

EO=SPECIFIC ENERGY AT INITIAL RADIUS RO 

M=S/C MASS (KG) 

MO=INITIAL S/C MASS (KG) 

ME=ARBITRARY MASS FOR NONDIMENSIONALIZING (20000KG) 
MBAR=NONDIMENSIONALIZED MASS 

MF=MASS OF FUEL BURNED IN TIME INCREMENT 

MFK=MASS OF FUEL BURNED PER TIME INCREMENT KEPLERIAN 
MFT=MASS OF FUEL BURNED TOTAL 

MFTK=MASS OF FUEL BURNED TOTAL WITH KEPLERIAN MOTION 
GAMMAR=FLIGHT PATH ANGLE (RAD) 

GAMMAD=FLIGHT PATH ANGLE (DEG) 

TH=THRUST (N) 

THK=KEPLERIAN THRUST 

THMAX=BLOWDOWN (MAX) THRUST (N) 
THBAR=NONDIMENSIONALIZED THRUST =(TH*TAU**2) / (MR) 
THBM=NONDIMENSIONALIZED BLOWDOWN (MAX) THRUST 
ALPHAR=THRUST ANGLE (RAD) 

B=BALLISTIC COEFFICIENT (M/(CD*S) ) 
BBAR=NONDIMENSIONALIZED B =B/(RHOO*RE) 
RHO=CALCULATED ATMOSPHERIC DENSITY 

RHOO=REFERENCE ATMOSPHERIC DENSITY 

SPI=SPECIFIC IMPULSE 

V=VELOCITY 

VBAR=NONDIMENSIONALIZED VELOCITY 

G=GRAVITATIONAL ACCELERATION 

SMA=SEMI-MAJOR AXIS 

ECC=ECCENTRICITY 

T=TIME 

TINC=INCREMENT OF TIME(STEP SIZE) 
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C 
C 


TF=END TIME 
PTI=PRINT TIME INTERVAL 


CONSTANTS 


RE=RADIUS OF EARTH 
MU=EARTHS GRAVITAIONAL CONSTANT 
TAU=ORBITAL PERIOD AT RE 


START PROGRAM 
PROGRAM ORBMAINT 


VARIABLE DECLARATION 
IMPLICIT REAL*8(A-H,M-Z) 
DIMENSION XBAR(4),XBDOT(4) 


CONSTANT DEFINITIONS 
PI=DATAN(1.0D+00) *4.0D+00 
MU=3.98601208133D+14 
G=9.806D+0 

RE=6.3782D+6 
TAU=2.*PI*RE**1.5/MU**0.5 
J=4 

ME=20000.D+0 


INITIALIZE FUEL USAGE TOTALS 


MFT=0.0D+0 
MFTK=0.0D+0 


CCCCCCCCCCCECECCCCEC ECC CE CECECEC CECECE CCC ECE CCE EGE CE CEC ECE CEC CeCe 


c 
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MAIN PROGRAM 


OPEN (10, FILE=‘ inp’, STATUS=‘OLD’ ) 

OPEN (20, FILE=‘outp’ , STATUS=’ NEW’ ) 
OPEN (30, FILE=‘orbp’ , STATUS=‘NEW’ ) 
OPEN (40, FILE=‘ ratp’ , STATUS=’ NEW’ ) 


READ INPUTS AND INITIALIZE PARAMETERS 
READ (10,1)RO,MO, THMAX, TB, TFB, TINCB, PTIB, BBAR, BETA, SPI 
FORMAT (10(/,21X,D13.7)) 


M=MO 
MBAR=M/ME 
TAURO=2.*PI*RO**1.5/MU**0.5 
TINC=TINCB*TAU 
THBM=THMAX* TAU* TAU / (ME*RE) 


READ PARAMETERS THRUST ANGLE, DESIRED BAND,THRUST FACTOR 


PRINT*, ‘ENTER ALPHA’ 
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READ*, ALPHA 
MGEHAR=ALPHA*PI/180.0D+0 


PRive*,  SNIER BAND i) 
READ* , BAND 
BAND=BAND* 1000. 0D+0 


PRINT A) “ENTERGIAPAC © 
READ*~, THEAC 
THRPAC-THRAC*].0D+0 


INITIALIZE COUNTERS USED IN PROGRAM 


INDEX=0 
KOUNT=1 
TRAD=0.0OD+0 


INITIALIZE NONDIMENSIONALIZED STATE VARIABLES 


XBAR(1)=RO/RE 

XBAR(2)=OD+0 

XBAR(3)=OD+0 
XBAR(4)=2.0D+0O*PI*TAU/TAURO 


COMPUTE INITIAL ORBITS VELOCITY AND ENERGY 


VO=(MU/RO)**0.5 
EO=((VO*VO) /2.-MU/RO) 


COMPUTE RADIUS, VELOCITY,ENERGY FOR DESIRED BAND MIN AND MAX 


RMIN=RO-BAND/2.0 
VMIN=(MU/RMIN) **0.5 
EMIN=( (VMIN*VMIN) /2.-MU/RMIN) 


RMAX=RO+BAND /2.0 
VMAX=(MU/RMAX)**0.5 
EMAX=( (VMAX*VMAX) /2.0-MU/RMAX) 


SET UP OUTPUT FILE HEADERS 
IBAND=INT (BAND) 


ITHMAX=INT(THMAX) 
IALPHA=INT (ALPHA) 


WRITE(20,*)‘!!BAND=',IBAND,’ THMAX=',ITHMAX,’ ALPHA=', IALPHA 
WRITE(20,*)‘!!EMIN=',EMIN,’ | EMAX=‘,EMAX 
WRITE(30,*)‘!!BAND=',IBAND,’ THMAX=',ITHMAX,’ ALPHA=‘,IALPHA 
WRITE(40,*)‘!!BAND=’,IBAND,’ THMAX=',ITHMAX,’ ALPHA=’,IALPHA 
WRITE(20,*)’!! ORBITS RADIUS VELOCITY MFB MFBK 

* ANGM ENERGY ’ 

WRITE(20,*)‘!! (KM) (KM/SEC) (KG) (Ke) 
WRITE(30,*)‘!! ORBITS SMA EGC APOGEE PERIGEE 
x PERIOD’ 

WRITE(40,*)‘’!! TIME ORBITS DRAG THRUST MASS 

* GAMMA ‘ 
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CALCULATIONS 


CALL DRAG TO SET INITIAL VALUE FOR FORCED KEPLERIAN MOTION 
USED TO COMPUTE KEPLERIAN THRUST 


CALL DRAG(BBAR, RE, XBAR, BETA, MBAR, VBAR,D) 
DK=D 
TK=DK*MBAR*ME*RE/ (TAU*TAU) 


CALL SUBROUTINE DRAG TO COMPUTE PRESENT DRAG 
CALL DRAG(BBAR, RE, ABAR, BETA, MBAR, VBAR,D) 


CALCULATE E TO SEND TO THRUST SUBROUTINE TO DETERMINE 
IF THRUSTER SHOULD BE ON OR OFF 


E=(RE*VBAR/TAU) **2.0/2.0D+0-MU/ (XBAR(1)*RE) 


CALL SUBROUTINE THRUST TO DETERMINE IF THRUSTERS ARE CN@ee 
OFF AND THE SET THE VALUE OF THE THRUST TERM 


CALL THRUST(EO,XBAR, THBAR, THFAC, RMIN, EMAX,E,ETH, THBM) 
CALL SUBROUTINE EQN TO UPDATE ORBITAL EQUATIONS 
CALL EQN(D, XBAR, XBDOT, VBAR, THBAR, ALPHAR, MBAR, PI) 


CALL SUBROUTINE RK4 TO USE RUNGE-KUTTA FOR INTEGRATION 
OVER TIME 


CALL RK4(TB, XBAR, XBDOT,J,TINCB, INDEX) 
CHECK TO SEE IF RUNGE-KUTTA COMPLETE 


TP (INDEX, .NE- ©) Co ToOMeG 


RUNGE=KUTTA DONE FOR THIS TITERATION OR. ORS tT 
UPDATE VARIABLES FOR CALCULATION OF ORBITAL 
PARAMETERS AND FOR OUTPUT 


TH=THBAR*ME*RE/ (TAU*TAU) 
R=XBAR(1)*RE 
V=RE*VBAR/TAU 

MF=TH*TINC/ (SPI*G) 

M=M=-MF 

MFT=MFT+MF 
GAMMAR=ATAN (XBAR(2) /(XBAR(1)*XBAR(4) ) ) 
GAMMAD=GAMMAR*180.0D+0/PI 
ANGM=R*V*COS (GAMMAR) 
ENERGY=((V*V) /2.0)-(MU/R) 
TRAD=TRAD+TINCB 
OANG=TRAD*2.*PI 


COMPUTE KEPLERIAN FUEL VALUES 


MFK=TK*TINC/ (SPI*G) 
MF TK=MFTK+MFK 
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C CALL SUBROUTINE ORBPAR TO COMPUTE ORBITAL PARAMETERS 
Cc 
CALL ORBPAR(ENERGY,ANGM,R,MU,PI,ECC,SMA,APOGEE, PERIGE, PERIOD) 
S COMPUTE NEW MBAR AND PUT R,V IN OUTPUT FORM 
ic 
MBAR=M/ME 
R=R/1000.0 
V=V/1000.0 
Cc 
Cc DATA OUT 
Cc 
C CHECK TO SEE IF TIME TO PRINT OUTPUT 
Cc 
IF(KOUNT .LT. DNINT(PTIB/TINCB)) GO TO 200 
Cc 
Cc PRINT OUTPUTS TO OUTPUT FILES OUT, RAT, ORB 
Cc 
WRITE(20,2)TB,R,V,MFT,MFTK, ANGM,E 
p BORMEGW hCG ho 4) 20 HONe oan 490%, FS 22 02% ,.FS ee 2k, F1i2.0, 
PO aE. 0} 
C 
WRITE (30,3)TB, SMA, ECC, APOGEE, PERIGE, PERIOD 
3 PORMAm( Ging so. 1 Wi@ gear deg OX, F10,3,1X,Fi0.3,2%,Fe.4) 
Cc 
WRITE(40,4)T,TB,D, TH,M,GAMMAD 
4 MORMAT( 26840, 1% h602,9X%) 612.9, 3X,6500,5xX,F9.3,1%,F10.4) 
Cc 
Cc RESET COUNTER 
Ec 
KOUNT=0 
Cc 
Cc CHECK TO SEE IF ORBIT COMPLETE, RESET ORBITAL ANGLE 
C 
IF (TRAD .GE. 1.0D+0) TRAD=0.0D+0 
Cc 
c UPDATE COUNTER 
200 KOUNT=KOUNT+1 
c 
Cc CHECK FOR END OF PROGRAM 
Cc 
Pio eins CO. TOM O0 
C 
END 
S 
C 
ie ec eCCCCCCCECeCCeC CCC CeCCCCECcceccceceeccccccceccccceccecccceccecccc 
Cc 
c THIS SUBROUTINE CALCULATES DRAG FOR USE IN THE EQUATIONS 
C OF MOTION. 
Cc 


SUBROUTINE DRAG (BBAR, RE, XBAR, BETA, MBAR, VBAR, D) 
IMPLICIT REAL*8(A-H,M-Z) 


C 
DIMENSION XBAR(4) 
e 
Cc SET REFERENCE ALTITUDE FOR ATMOSPHERIC DENSITY CALCULATIONS 
Cc 
RREF=6.638145D+6 
Cc 
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C COMPUTE PRESENT NONDIMENSIONALIZED VELOCITY 
VBAR=( (XBAR(2) *XBAR(2))+(XBAR(1)*XBAR(4))**2)**0.5 


2. 
C COMPUTE EXPONENTIAL TERM IN DRAG EQN 
C 
E=EXP (-BETA*RE* (XBAR(1)-—(RREF/RE) ) ) 
C 
c COMPUTE DRAG TERM 
e 
D=E*VBAR*VBAR/ (2.0*MBAR*BBAR) 
e 
C 
RETURN 
END 
Ce: 
@: 
CCCCCCCCCECECCEGCCCECECECECCEECECEECECCCCOCcCCeeECeeeereeeeee 
eC 
@: THIS SUBROUTINE UPDATES THE EQUATIONS 
G OF MOTION 
C 


SUBROUTINE EQN(D, XBAR, XBDOT, VBAR, THBAR, ALPHAR, MBAR, PT) 
IMPLICIT REAL*8(A-H,M=Z) 


e 
DIMENSION XBAR(4),XBDOT(4) 
e 
Cc 
A=XBAR(1)*XBAR(4)*XBAR(4) 
B=4.0*PI*PI/(XBAR(1)*XBAR(1)) 
C=D* (XBAR(2) /VBAR) 
E=THBAR*SIN(ALPHAR) /MBAR 
Cc 
Cc 
F=2.0*XBAR(4)*XBAR(2)/XBAR(1) 
G=(D/XBAR(1))*(XBAR(1)*XBAR(4) /VBAR) 
H=THBAR*COS (ALPHAR) / (MBAR*XBAR(1) ) 
Cc 
6 
XBDOT (1)=XBAR(2) 
XBDOT (2)=A-B-C+E 
XBDOT (3)=XBAR(4) 
XBDOT (4)=-F-G+H 
@ 
e 
RETURN 
END 
Cc 
é 


ccoceccecceccccccececcececececccccceeesecacees ceeerececceeeecces 
Ee 
C THIS SUBROUTINE CALCULATES ORBITAL PARAMETERS 
C 
SUBROUTINE ORBPAR(ENERGY,ANGM,R,MU,PI,ECC,SMA,APOGEE, PERIGE, 
* PERIOD) 
IMPLICIT REAL*8(A-H,M-Z) 


PROBLEM COULD ARISE IN COMPUTING ECCENTRICITY IF VERY 
SMALL- COULD TRY TO TAKE SQRT OF NEG NUMBER. TO 
PREVENT, IF ECCENTRICITY IS LESS THAN 1E-6, SET 
EQUAL TO 0. 


OO Aa:0°80:Q 
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PROB=(1.12.*ENERGY *ANGM*ANGM / (MU*MU) ) 
Poo DobS (PROB) sol. 1. D=12). THEN 
ECC=0.0D+0 
ELSE 
ECC=PROB**0.5 
oONDIF 


SMA=(-MU/(2.*ENERGY))/1000.0 


SMAM=SMA*1000.0 

APOGEE=SMA* (1.0+ECC) 

PERIGE=SMA* (1.0-ECC) 

PERIOD=2.0*PI* (SMAM**1.5)/(MU**0.5)/60.0 


RETURN 
END 


PREPS CCCEe CeCe CCCCEC CCC CC CECCCCC EC CECE CECCCECECE CCCCCCCCECCCCCECCCE 


THIS SUBROUTINE DOES A FOURTH LEVEL RUNGE-KUTTA 
INTEGRATION. 


SUBROUTINE RK4(TB,XBAR, XBDOT, J, TINCB, INDEX) 
PMP olCil REAL*S(AqH, M=Z) 
INTEGER INDEX,1 


DIMENSION XBAR(4),XBDOT(4),SAVED(4),SAVEX(4) 


INDEX=INDEX+1 
GO TO (1,2,3,4), INDEX 


DO 10 I=1,J 

SAVEX (I) =XBAR(TI) 

SAVED (I) =XBDOT (I) 
XBAR(I)=SAVEX(I)+0.5D+0*TINCB*XBDOT(TI) 
TB=TB+0.5D+0*TINCB 

RETURN 


DO 20 I=1,J 

SAVED (I) =SAVED(I)+2.D+0*XBDOT (TI) 
XBAR(I)=SAVEX(1I)+0.5D+0*TINCB*XBDOT (I) 
RETURN 


DO 30 I=1,J 

SAVED (I) =SAVED (1) +2.D+0*XBDOT (I) 
XBAR (I) =SAVEX(1)+TINCB*XBDOT (I) 
TB=TB+0.5D+0*TINCB 

RETURN 


DO 40 I=1,J 

XBAR (I) =SAVEX(1I)+TINCB/6.D+0* (SAVED (I)+XBDOT(I) ) 
INDEX=0 

RETURN 
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Cc 
C 


END 


CCCCCCCCCCCCCCCCCCECCOCCCCCCECECCC CCCCCECCC CCC CCC CEC CEC CEC CECE See 


THIS SUBROUTINE DETERMINES IF THRUSTERS SHOULY SE Cie. 
OFF TO TRY TO STAY WITHIN THE DESTREDS aA 


SUBROUTINE THRUST(EO, BAR, THEAR, THYAG, RMIN Man], 2 


*THBM) 


TMPLICIT REAL*8(A-H,M-Z) 
DIMENSION XBAR(4),XBDOT(4) 
RE=6.3782D+6 


COMPUTE NONDIMENSIONALIZED POINT FOR THRUSTERS TO TURN ON 
BASED ON RADIUS. 


RTH=RMIN/RE 
ETH=EMAX~— (EMAX-EMIN) * THFAC 
IF THRUSTERS ARE ON GO TO SEE IF TIME TO TURN OFF 
IF (THBAR .EQ. THBM) GO TO 100 
CHECK IF RADIUS BELOW BAND 
IF (XBAR(1) .LE. RTH) THEN 
CHECK IF ENERGY BELOW INITIAL ENERGY 
IF (E .LE. EQ) THEN 
TURN THRUSTERS ON 


THBAR=THBM 
ELSE 


OTHERWISE LEAVE OFF 


THBAR=0D+0 
ENDIF 


RADIUS NOT BELOW BAND, LEAVE OFF THRUSTERS 
ELSE 
THBAR=0D+0 
ENDIF 
IF THRUSTERS ARE ON, SEE IF TIME TO TURN OFF 
IF (E .GE. EQ) THBAR=OD+0 


RETURN 
END 


ccceccccceccececccecccccccccecceccececccececcccecccccececcecceccccccccccccece 
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B. SAMPLE INPUT FILE 


C. SAMPLE OUTPUT FILE 


''BAND = 
'EMIN= 


1.00 
2.00 
3.00 
4.00 
5.00 
6.00 
7.00 
8.00 
9.00 
10.00 
EOO 
12.00 
700 
14.00 
15.00 
16.00 
17.00 
18.00 
19.00 
20.00 
21.00 
e200 
23.00 
24.00 


25000 THMAX= 
-30080181.48670808 
'! ORBITS RADIUS 


(KM) 
6636.385 
6634.580 
6632.737 
6630.855 
6628.930 
6626.95 1 
6634.414 
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